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ABSTRACT 

We examine various physical processes associated with the formation of globular clusters by using the 
three-dimensional Smoothed Particle Hydrodynamics (SPH) code. Our code includes radiative cooling of 
gases, star formation, energy feedback from stars including stellar winds and supernovae, and chemical 
enrichment by stars. We assume that, in the collapsing galaxy, isothermal cold clouds form through 
thermal condensations and become proto-globular clouds. We calculate the size of proto-globular clouds 
by solving the linearized equations for perturbation. We compute the evolution of the inner region of 
the proto-cloud with our SPH code for various initial radius and initial composition of gases. When the 
initial gases contain no heavy elements, the evolution of proto-clouds sensitively depends on the initial 
radius. For a smaller initial radius, the initial star burst is so intense that the subsequent star formation 
occurs in the central regions to form a dense star cluster as massive as the globular cluster. When the 
initial gases contain some heavy elements, the metallicity of gases affects the evolution and the final 
stellar mass. If the initial radius of the proto-globular clouds was relatively large, the formation of a star 
cluster as massive as the globular clusters requires the initial metallicity as high as [Fe/H] > — 2. The 
self-enrichment of heavy elements in the star cluster does not occur in all cases. 

Subject headings: globular clusters: general — galaxies:star clusters — stars: formation — 
hydrodynamics 



L INTRODUCTION 

The globular clusters belong to the oldest populations 
m o ur galaxy. For the general rcvic wa on globular oluatcra, 
sec |Moylan, fc Hoggio (1997D and [Harris (199l| ). Their- 



lar clusters form from large-scale unorganized mo- 
tion of interstellar gas as in the Magellanic Clouds 
which are now producin g young clusters (Kumai 



Basu, & Fujimoto 1993). Another example is that 



forraation ( Elmergreen et al. 1999| ) is closely related to the 
formation process of our galaxy. The formation of a globu- 
lar cluster may take place in two stages: (1) the formation 
of a proto-globular cloud (PGC) and (2) the formation of 
a star cluster from the PGC. There are three scenarios for 
the formation of a PGC, denoted as primary, secondary, 

where the PGC 
or after the 



very young globular clusters are observed in nearb y 
galaxy NGC 1705 and 1569 (|Ho fc Filippenko 19961). 



and tertiary model (Fall & Rees 1987), 
forms in different stages, i.e., before, during, 
collapse of the galaxy, respectively. 



The primary model was suggested by Peebles & 



Dicke 

of the recom 



(1968) 



They showed that the Jeans mass 
Dination stage of the universe is com- 
parable to the observed masses of globular clusters 
and thus the PGC can form due to the gravitational 
instability. Globular clusters may be debris of these 
objects. The most serious problem with this primary 
model is that there have been few intergalactic glob- 
ular clusters discovered. Almost all globular clusters 
ever discovered exist in galaxies. 

In the secondary model, on which we concentrate in 
the present paper, the PGC fo rms due to thermal in- 
stabilities ( Fall fc Rees 1985| ). Our detailed investi- 



Once the PGC forms, there might be many ways to form 
a globular cluster (GC). Murray fc Lin (1993) summarized 
the scenarios of the formation of GCs from PGCs as fol- 
lows. The PGC can be divided into two types depending 
on their masses. The cloud, whose mass exceeds the Jeans 
mass, is gravitationally unstable, thus spontaneously col- 
lapsing to form stars. The cloud, whose mass is smaller 
than the Jeans mass, is stable until some instabilities are 
introduced. A cloud-cloud collision or cloud-disk collision 
can trigger such instability. When such collisions occur, 
the cloud would become thermally unstable owing to the 
efficient cooling. This cooling would lead to the forma- 
tion of a very dense region, thus inducing a burst of star 
formation. 

In the present paper, we examine quantitatively the 
above still qualitative scenario of the globular cluster for- 
mation. We use our Smoothed Particle Hydrodynam ics 
(SPH) code ([Lucy 1977|; [Gingold fc Monaghan 1977|) to 



gation and calculation are presented in the following 
sections. 



For the tertiary model. 



one example is that globu- 
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simulate the formation of a globular cluster from a PGC. 
Our code includes the following physical processes: radia- 
tive cooling, star formation, energy feedback from stars in- 
cluding stellar winds and supernovae, and chemical enrich- 
ment from stars. The SPH method which includes star for- 
mation processes like ours has been applied to many astro- 
physical problems. Such problems include the formation of 
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isolated galaxies ( Katz 1992 : Stcinmetz fc Miiller 1994 ), 
the evolution of galaxies ( Friedli fc Benz 1995|), a nd the 
cosm qlogical simulations ( Navarro fc White 1993 ). How- 



ever, 3ur study is the first attempt to apply this method 
to the globular cl uster formation (prelim i nary results have 
been presented in Nakasato et al. (1996|); Nakasato, Morij 



& No jnoto (1999^ 



Among the two triggers of the instability, collapse and 
collision, we concentrate on the collapse case in the present 
paper. The collision case for a wide range of parameters 
(masses of the two clouds, a relative velocity between the 
clouds etc.) will be discussed elsewhere. This paper con- 
sists of the following sections. In section 2 our SPH method 
is described. In section 3, we present our linear analysis 
for the thermal instabilities to form a PGC. In section 4, 
we present the results of the evolution of the PGC with 
two different initial compositions of gases. Sections 5 and 
6 are devoted to discussions and concluding remarks, re- 
spectively. 

2. METHOD 

To simulate the formation of a star cluster from gases, 
we use our GRAPE-SPH code using Rc mote-GRAPE li- 
brary ( [Nakasato, Mori, fc Nomoto 19971 ). The GRAPE is 
a special purpose computer for the self gravity calcu lations 
in general N-body system ( Sugimoto et al. 1990| ). Our 
GRAPE-SPH code includes various physical processes, 
e.g., radiative cooling, star formation, and feedback from 
formed stars. There are many different implementations 
of the SPH method (a recent co mparison of various SP H 
implementations is presented in Thacker et al. (1998)). 



In the present paper, we describe an essential point of our 



SPH (Tude. The SPH fuiniidaLiuu LhaL we use is Lhe same 
as th e described in 



Navarro fc White (19 93 ). Wc use 



the smoothing length that can vary spatially and evolve 
with time, and integrate equations of motion with a sec- 
ond order Run ge-Kutta method as described in Navarro & 



White (1993). Details of the imple mentation of our code 
are described in Mori et al. (2000). In the following sub- 



sections, we describe the physical processes (cooling, star 
formation, and feedback) in some detail. 

2.1. Radiative cooling 

The radiative cooling rates depend on the temperature 
and ionization states of the gas. Also, the chemical com- 
position of the gas affects the cooling rate drastically. We 
perform SPH simulations for the following three cases (A, 
B, C) of the gas with different physical state and chemical 
composition. 

Case A: We assume that the chemical composition of 
the gas is primordial with no heavy elements and the gas 
is in ionization equilibrium. In this case, our treatment 
of radiative cooling is essentially the same as adopted by 
Katz, Weinberg, fc Hernquist (1995| ); they computed the 



cooling rate using the two-body processes of H and He, and 
free-free emissions. The cooling function (A(T)) is shown 
in Figure |l| with the solid line. In this case, the cooling 
rate decreases very rapidly as the temperature T decreases 
below 2 X 10^ K so that the gas would not radiatively cool 
below T - 10" K. 

Case B: We assume that the gas includes some heavy 
elements and is in ionization equilibrium. In this case, 
we use the cooling function with different chemical com- 



position that is computed by MAPPINGS HI software 
by R.S. Sutherland (MAPPINGS HI is the successor of 
MAPP INGS II that is described in [Sutherland fc Dopita 
(1993| )). We compute the cooling function of the ioniza- 
tion equilibrium gas for [Fe/H] = — 5.0 - 0.0 with the so- 
lar ab undance ratio (see Table 4 of Sutherland fc Dopita[ 
(1993| )) and present the results in Figure ^with the dashed 
lines (each line corresponds to [Fe/H] = — 1, —2, —3 respec- 
tively from the top to the bottom). Existence of heavy ele- 
ments significantly enhances the cooling rates. At T < IC* 
K, the cooling due to the forbidden line emission of car- 
bon and oxygen is efficient. For [Fe/H] = — 1, the cooling 
rates around T lO^K are 100 times larger than the cool- 
ing rates of primordial hydrogen and helium gas. These 
differences would make the evolution of the gas very dif- 
ferent. 

Case C: We concern about the non-equilibrium cooling. 
If the gas cools from high temperatures, ionization equi- 
librium is not realized. Figure | shows the ratio between 
the recombination time (trccom) of hydrogen and the cool- 
ing time (tcooi) of primordial gas in ionization equilibrium; 

^recom/^cool- The trecom 

and tcooi are defined as 



1.0 



and 



cool 



aiT)fe{T) 



3kbr 
2A(T) ' 



(1) 



(2) 



where the a is the recombination coefficient of hydrogen 
and electron, /e is the fraction of free electrons, a nd kb is 
the B olzman constant. For a, we use the value of 3pitzei 
(197^ ). Clearly, the cooling time is much shorter than the 
recombination time for T ~ 2 x 10^ K, so that ionization 
equilibrium is not realized when the gas cools from high 
temperatures. In the non-equilibrium case, the existence 
of electrons and ionized hydrogen at T < 10* K makes 
it possible to form H2 molecules through the creation of 
intermediaries H^ and HJ as: 



and 



H + e ^ H + photon 
H + H- ^ H2 + e 



H + H+ ^ H+ + photon 



H 



H 



Hs + ff 



(3) 



(4) 



( Shapiro fc Kang 1987 ). These H2 molecules cause fur- 
ther cooling down to T ~ 10^ K owing to the rotational- 
vibrational line excitation. 

To include the effect of such molecular cooling, we have 
to solve rate equations that determine the ionization states 
of H and He atoms, and the formation and destruction of 
H2 molecules. In our case, the typical time step in solving 
the rate equations is shorter than the dynamical time step, 
which is determined mainly by the Courant condition. In- 
tegrating all equations with a shorter timestep than the 
dynamical time is very costed work. So we divide the dy- 
namical timestep with dynamical variables (density etc.) 
being constant in solving the energy an d rate equations, 
which is similar to the method adopted in ^hapiro fc Kangj 
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(1987). The rate-coefScie nts are also the same as used in 
Shapho fc Kang (1987| ). Included species are H*^, H+, 
He°, He+, He+ + , H", H+, H2, and e, where is the 
exited hydrogen molecule. In our SPH code, we solve the 
rate equations for 10 species in each SPH particle with a 
reasonable computing time. Solving the rate equations for 
the gas including heavy elements (over 200 species) is not 
feasible with a current resource so that we only concern 
the hydrogen and helium plasma in the present paper. 

The cooling function for the non-equilibrium case is pre- 
sented in Figure H with the dotted line. In computing these 
cooling rates, we follow the isobaric temperature evolution 
of the fully ionized gas. Initially, the temperature and hy- 
drogen number density of the gas are 10^ K and 0.01 cm"'^, 
respectively. We assume optically thin plasmas so that the 
gas cools rapidly. Around T ~ 2 x 10''' K, H2 molecules 
begin to form. At T < 10^ K, the cooling rate due to H2 
molecules is comparable to the cooling rate for [Fc/H] = 
-1 gas. 

Finally, we summarize the treatment of the cooling rates 
in our SPH code. We can perform SPH simulations for 
these three cases A, B and C. In the present paper, we 
will study the evolution of a PGC with two different ini- 
tial chemical compositions, e.g., a metal- free gas cloud and 
a metal-rich gas cloud. In the former case, we solve the en- 
ergy and rate equations simultaneously for each SPH par- 
ticle (case B) if the gas temperature is lower than 3 x 10"* K 
and the gas particle is not in the heating phase (see section 
2.3). If these conditions are met, we use the pre-computed 
cooling table for the ionization equilibrium case (case A). 
In the latter case, we use the pre-computed cooling table 
to solve the energy equations (case C). 

2.2. Star formation 

Our treatmen t of star formation is the same as adopted 
by lKatz (199^ ). Hereafter, "STAR" means "star parti- 



cle" , which is not an individual star but an association of 
many stars. A STAR forms in the region where the flow is 
converging, cooling, and Jeans unstable. These conditions 
are expressed as 

(V • v), < 0, (5) 



^cool ^ ^d; 

td < ^sound ■ 



(6) 
(7) 



Here icooi, , and isound are the cooling time, dynamical 
time, and sound crossing time, respectively, and expressed 
as 



^cool 



1 



t 



sound 



(8) 
(9) 
(10) 



where /x is the mean molecular weight, u is the specific 
thermal energy, and Cs is the local sound speed. A STAR 
forms in the region where these all three conditions are 
satisfied. The star formation rate is given as 



Dp* 
Dt 



''Star form 



(11) 



where p» is the density of a star, and istarform = td/C with 
a star formation parameter C. In the present paper, we 
adopt C = 1.0, i.e., a STAR forms in a dynamical time. 
The last term of eq uation p] is similar to the Schmidt's 
law (Schmidt 1959). Integrating equation (0) over one 



time step 5t and making calculations with the equations 
for SPH, we obtain the mass of a newly formed star in St 
as 



^star 



1 



1 + 0.5 



St 



t=tarfo 



■k\ Pi 



(12) 



The newly formed STAR is then treated as a collision-less 
particle. We introduce the minimum allowed mass (~ 10 
Mq) for the newly formed STAR to prevent the unphysi- 
cal effects in the equation of motion and the treatment of 
the feedback processes. 

We note the star formation recipes in our SPH code. In 
the regions with increasing densities, first two conditions 
for the star formation (Eq. ^ and ^ are almost satisfied. 
Thus whether STAR forms in some regions is determined 
mostly by the Jeans criterion (Eq. The critical density 
for the star formation in our treatment is obtained as 



PJcans 



\4:7rGpm J m~ \/3' 



(13) 



Here, it is assumed that the neighbor radius is expressed 

as 



1/3 



(14) 



where /3 is determined experimentally, because hi is deter- 
mined in order to make the number of neighbor particles 
almost constant in some range (30 - 80 in our codes). Since 
the typical value of /3 is 1.0 - 1.1, (3^ ranges 1.0 - 2.0. Then 
if rrii is constant, the critical density is determined almost 
only by the temperature. The typical calculations in the 
present paper use 5000 gas particles for a 10^ Mq gas 



sphere. Thus the initial mass per particle (m^ in Eq. 13) 
is ~ 200 Mq. For T = lO"* K and 10^ K, the critical 
density for star formation is pjoans ^ 10^^^ and 10^^^ g 
cm^^, respectively. With the star formation recipes used 
in our SPH code, the STAR forms only in the very high 
density region if T - lO'' K. There are, however, the max- 
imum density (pmax) that numerically reached in the SPH 
method; it is estimated 



Nnrrii 



(15) 



where is the number of neighbor particles and e is the 
gravitational softening length. In the present calculations, 
Pmax ^ 1-8 X 10"^'^ g cm~'^. These arguments ensure that 
with the star formation recipes used in our SPH code and 
the initial conditions of the present calculations, the STAR 
forms in the region where the temperature is as low as 10^ 
K. 
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2.3. Feedback 

The formed stars eject gases and heavy elements in stel- 
lar winds and Type II supernova explosions and heats up, 
accelerate, and enrich a circumstellar and an interstellar 
medium. High energy explosions like supernova produce 
high temperature and low density regions in interstellar 
medium. In the SPH method, the numerical accuracy for 
high d ensity regions is much better than mesh baaed mcth - 



Mq are assumed. In Eq. (17) and (18), M,„a (= 50.0 
Mq) and {— 8.0 M©) are the upper and lower hmit 
masses of the stars that explode as Type II supernovae. 

For the supernova energy, we assume that esNii = 10^^ 
erg. For the stellar wind, esw is estimated to be 0.2x10^^ 
erg for solar metallic ity stars from the observational data 
of OB associations ( Abbot 1982| ). The Chemical abun- 



dance of a massive star sig nificantly affects esw (Leitherer 



ods b mt the accuracy for low density rogione are p oorer. A 



typical resolutio n of usual SPH simulations, e.g., Navarro 



& White (1993), including a star formation process (100 



1000 pc) is larger than a typical size of supernova rem 
nants (< 100 pc). Thus, it is difficult to properly include 
the energy, momentum, and mass release from stars in the 
SPH method because of the nature of the SPH method 
and the poor resolutions in current computing resources. 
We must use some approximations to mimic real feedback 
processes in the SPH method 



One of the method has been proposed by Navarro & 
White (1993). In their method, the energy produced by 



a supernova explosion is distributed to neighbor gas par- 
ticles of each STAR mostly as a thermal energy and the 
rest is distributed as a velocity perturbation to the gas 
particles; the fraction of the energy in a kinetic form is a 
free p arameter (we note that Leitherer, Robert, fc Drissen] 
(1992| ) presented the population synthesis model of stellar 



feedback processes). In the present paper, we distribute 
the energy to neighbor particles in a pure thermal form as 
a zero-th order approximation, because the size and time 
scale of our model a re much smaller than thos e of a galaxy 
formation model of Navarro & White (1993). 



2.3.1. Energy ejection 
The energy ejection rate per STAR is given as 



'eject 



esw^sw + esNii^SNii, 



(16) 



where esw is the total energy ejected by stellar winds dur- 
ing the stellar lifetime and esNii are the energy ejected by 
one Type II supernova explosion. The i?sw is the number 
of stars per unit time expelling their envelopes at the cur- 
rent epoch and i?sNii is the rate of Type II supernovae. 
We define the i?sw and i?sNii as follows 



^sw = 



-RsNii 



A/up 



^{m)dm 



t{M„ 



^{m)dm 



M„ 



-(Af„ 



T(M^a) ' 



(17) 



(18) 



where $(m) is the initial mass function (IMF), namely 
<^{rn)dra gives the number of stars in the mass r ange of 
(m, m + dm) and the r (m) is the stellar lifetime (David. 
Forman, fc Jones 199C| ). In the present study, we assume 
the power law type IMF as 



$(r7i) — Am 



'2.35 



(19) 



where the A is the constant. For the upper and lower 
limit masses in Eq. (19), M^^ = 120 Mq and Mio = 0.05 



Robert, & Drissen 1992), so that we use metallicity de- 
pendent esw as esw = 0.2esNii(2'/ZQ)°-^, 
mass fraction of metal in the STAR. 



where Z is the 



2.3.2. Mass ejection 

In our code, the mass ejection due to stellar winds is 
combined with the mass ejection due to Type II super- 
nova. Thus, the mass ejection rate per STAR is written 
as 

Afcjcct = "^SNII^SNII, (20) 

where msNii is the average mass ejected by stellar winds 
and Type II supernovae defined as 



msNii 



TO$(m)(im 



M„ 



AI„- 



(21) 



<^{m)dm 



M„ 



Here ttins is the mass that is locked up in the neutron star 
and assumed to be 1.4 Mq. The fraction of heavy metal in 
Mcject is co mputed by the nucleosynthesis yield of Type I I 



supernovae (Tsujimoto et al. 1996; Nomoto et al. 1997) 



2.3.3. Summary 

We assume that the feedback phase continues for 
T(Mms) = 4.3 X 10'' yr from the formation of each STAR 
and is divided into two phases: a stellar wind phase and 
a supernova phase. The stellar wind phase continues for 



r(M„ 



5.4 X 10 yr, during which only the energy ejec- 



tion from STARs is included; the ejected mass is included 
in the supernova phase for simplicity. The supernova phase 
begins at t = T(M„ia) and ends &i t = t{M^s)- During 
the supernova phase, the energy ejection is given by Eq. 
(16). The mass ejection is the sum of the contributions by 
the stellar winds and Type II supernovae. 

The thermal energies, gases, and heavy elements from 
stellar winds and supernovae are smoothed over neighbor 
particles of the STAR within a neighbor radius of i?/ (feed- 
back radius). We treat i?/ as a parameter to meet the ob- 
servational constraints. Such neighbor particles are called 
"in heating phase" . When the gas particles are in heating 
phase, w e assume that the cooling is su ppressed as pro- 
posed in Mori, Yoshii, fc Nomoto (1999 ). This treatment 
produces the high temperature region around the STAR. 
Thus, the star formation is forbidden in the gas particles 
in heating phase. 

3. PROTO-GLOBULAR CLOUD FORMATION 

We first examine the radiative condensations, which oc- 
cur in a wide range of astrophysical ci rcumstances fro m so- 
lar prominence to interstellar clouds (Meerson 1996). Ra- 
diative condensations in optically thin plasma have been 
considered by many authors since the pioneering work by 
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Parker (1953 ) and Field (1965 ) . The scale length of grav- 
itational instability in a collapsing proto- galaxy is much 



larger than the radii of globular clusters (Lin & Murray 
1996^] As will be shown in the following sections, the scale 
length of radiative condensations is comparable to the radii 
of globular clusters. Thus, in a collapsing proto-galaxy, ra- 
diative c ondensations may be the only mechanism to form 
a PGC ( [Fall fc Rees 1985| ; |Lin fc Murray 1996| ). 



3.1. Radiative condensations in a collapsing proto-galaxy 

The characteristic equation for the growth rate of ra- 
diative condensations, n, is obtained from the linearized 
equations for perturbations as 



kx + - — +nc 



(22) 

where Cg is the sound speed, 7 is the ratio of the specific 
heats, k — 27r/A is the wavenumber of the perturbation, 
kp and kx are the wavenumber of sound waves whose fre- 
quencies are equal to the growth rate of isothermal and 
isochoric perturbation, respectively, and fcx is the invers e 
of the scale length of thermal conduction ( Field 1965|) . 
fcp, kx and k^ are expressed as 



kr = 



m(7 - i)Po Hn) 

ckbTo tJ? 

^(7 - 1) po dA 
ckb /i^ dT ' 



kK = 



m(7- 



(23) 



(24) 



(25) 



where Tq and po are the equilibrium temperature and den- 
sity, respectively, k^ is the Bolzman constant, and A(T) is 
the cooling function (Figure |^). We assume 7 = 5/3 and 
K = 5.6 X 10"^ T^-^ erg s"^ K"^ cm'^ Solving Eq. 
(|2^ ) as a cubic equation of n for different fc, we obtain the 
dispersion relation between n and fc. In applying to our 
galaxy we adopt Tn ~ 1.0 x 10^ K and po - 1.7 x IQ-^^ g 



cm ( Fall fc Rees 1985 ). With these values, the disper- 



sion relation has a peak at some fc . In the present paper, 
we assume that the scale for the maximum growth rate is 
typical scale of a PGC. To obtain the typical scale of a 
PGC for different Tq and po, equatio n (p^) is view ed as a 
quadratic in fc^ (see Section II (d) in [Field (1965| )). 



The results for different Tq and po are presented in Fig- 
ure ^. In applying to our galaxy, we obtain the scale length 
of ~ 600 pc. This scale length is consistent with the fol- 
lowing simple estimate. For the adopted cooling function, 
the wavelengths for kp and kx with Tq ~ 1.0 x 10^ K and 
po ~ 1.7 X 10~^^ g cm~'^ are respectively obtained as 

Xp > 10^ pc and < 1 pc. (26) 

This implies that the perturbation with a scale greater 
than lO'^ pc is dumped owing to the limit of the sound 
speed, and the perturbation with a scale smaller than 1 pc 
is also dumped by thermal conduction. Thus the typical 
scale of a PGC in our galaxy is estimated to be several 
hundreds pc and the mass of a PGC ranges from 10*" to 
10^ Mq. 



3.2. Density profile of a PGC 

The estimated scale of a PGC is larger than the present 
radius of globular clusters (10 - 100 pc). This implies that 
during radiative condensations, a PGC shrinks before star 
formation begins or star formation in a PGC occurs in 
the central region of the cloud. If we consider the metal- 
free PGC where the ionization equilibrium is achieved, the 
cooling time, tc, of the PGC is much shorter than the dy- 
namical time, td {tc/td ^ 0.25 for Tq ~ 1.0 x 10^ K and 
Po ~ 1.7 X 10~^^ g cm~^). After the short period of radia- 
tive condensations (~ tc), the PGC becomes a warm dense 
cloud (T ^ 10'^ K) surrounded by a hot thin gas with Tq 



(FaU fc Rees 1985) 



The qualitative discussion on the evol ution of such a 
metal -free PGC has been presented by Lin fc Murray 
1996| ). They argued that in the first stage of collaps- 



ing galaxy, the assumption of ionization equilibrium is not 
valid because of little radiative emission. If the ioniza- 
tion equilibrium is not achieved, the temperature of high 
density regions can be as low as ~ 10^ K due to hydrogen 
molecular cooling. In this case, the high density cold cloud 
with T — 10^ K is surrounded by a warm gas. The criti- 
cal mass of the isothermal sphere confined b y an external 



pressure depen ds on the cloud temperature (Ebert 1955 
[McCrea 1957 ). If the temperature of the clo ud is ~ 10^ K 
the cr itical mass is ~ 10^ Mq. According to Lin fc Murray 
|(1996| ), such a small cloud as M '--^ 10^ Mq collapses to 
make the first stars in the proto-galaxy. These first gener- 
ation stars are the source of radiative emission to maintain 
the ionization equilibrium. If such radiation continues long 
enough, a cold PGC will evolve almost isothermally with 
T ^ 10"* K and form the isothermal profile of p (x r~^. The 
core size of the PGC will decrease from the initial size of 
several hundreds pc to a size comparable to the observed 
globular cluster. According to these arguments, we can 
assume that the density structure of a PGC has a form of 
p oc r~^. 

For the initial conditions of the PGC, we use the King 
profile with fixed mass of M ~ 10^ Mq which nearly 
equals to the Jeans mass of the current conditions {po and 
To). We assume that this sphere represents the inner re- 
gions of the PGC. We examine three different cases for a 
initial radius of Ri — 150, 200, and 300 pc. The initial gas 
density profiles of the three cases are presented in Figure 
m For the smaller Ri , the initial concentration of the in- 
ner region of the PGC is higher. The temperature of the 
sphere is assumed to be T '-^ 10^ K. The properties of the 
inner region of the PGC, which are the initial conditions 
of our calculations, are summarized as follows: 

• The mass of the clouds is M 10^ Mq. 

• The radius of the clouds is Ri ^ 150 - 300 pc. 

• The clouds is an isothermal sphere with T ^ 10^ K. 

• Initially, the velocity of the cloud is zero. 

In the following calculations, the initial number of gas par- 
ticles is ^ 5000 so that the initial mass of the one gas par- 
ticle is ^ 200 Mq. The dependence of varying the initial 
number of gas particles is discussed in Section 4.3. In all 
cases, the gravitational softening length for all particles is 
set to be 1 pc and fixed during the calculations. 
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4. RESULTS 

First, we note on the feedback radius (i?/). We use the 
Stromgren radius (-Rst) of a typical OB star as Rf. The 
typical value is Rs t 10 - 100 pc, w here the density of the 
ISM is 1 cm^^ ( psterbrock 1974| ). The radius depends 
on the density of the ISM as i?st c>c n~^/^. In the central 
region of our initial models, n ranges from 100 cm~'^ to 
1000 cm^"^. Thus, we can estimate i?/ ~ 1 — 5 pc (using 
the largest value of Rst)- In the following two subsections, 
we choose Rf = 3 pc for all gas particles. The results 
for both cases are summarized in Ta ble h| . The results for 
different Rj are described in Section [Op 

4.1. Evolution of the metal-free PGC 

In this section, we describe the evolution of the metal- 
free PGC for case C. 

As presented in the previous sections, the metal-free 
PGC may evolve isothermally with T ~ 10'' K. The 
isothermal evolution of PGC is terminated when the ioniz- 
ing radiations from the first generation stars stop. With- 
out ionizing radiation, the PGC cools efficiently by II2 
molecules. Our calculations start from that time when 
the PGC begins to cool. We assume M = 10® A/q, and 
the initial temperature of T = 10* K and the King pro- 
file sphere. Since the size of the PGC at this stage is 
unknown, we examine three different cases for an initial 
radius of Ri = 150, 200, and 300 pc. 

In all three cases, the initial density of the central re- 
gion is high enough for efficient cooling to occur, so that 
the temperature of the central region decreases rapidly to 
~ 10^ K. Due to the high density and low temperature 
{r-^ 10^ K), a burst of star formation occurs in the cen- 
tral region in all cases. After the first star formation, the 
central region becomes the heating region and the temper- 
ature of the central region increases gradually. 

At t ~ 6 Myr, Type II supernovae start to occur and 
the temperature of the heating region increases rapidly to 
T ~ 10® K. The left panel of Figure ^ shows the evolution- 
ary changes in the central temperature for Ri — 150 pc. 
Such a high temperature region expands and the central 
density begins to decrease (see the right panel of Figure 

H) . The expansion of the central region leads to the for- 
mation of a shell structure as seen in the evolution of the 
gas density profile (Figure ^. The density profile at t = 
6 Myr, clearly shows the shell structure. After that time, 
the shell expands outwardly (see the right panel of Figure 

I) . 

The star formation after the shell formation occurs not 
in the central region but in the shell. Figure |^ show the 
change in the radius of the star forming region as a func- 
tion of time for Ri = 150 pc. We can see that the star 
forming region moves outward. At t = 10 Myr, where we 
stop the computation for Ri = 150 pc, the stellar mass 
reaches ~ 1.3 x 10^ AIq. This means that the star for- 
mation efficiency is 13 %. The bound stellar mass at 
t = 10 Myr is ~ 10^ Mq and the mass is comparable 
to the typical mass of globular clusters (~ 10^ Mq). At 
t = 10 Myr, the gas is almost removed from the central 
region. Figure |^ shows the projected particle positions of 
STAR particles (left panel) and the stellar density profile 
(right panel) at t = 10 Myr. The STAR particle shows 
elongated shape (bar like shape). Such a shape is caused 



by the radial orbit instability of the STARs formed at the 
large radius (Palmer & Papaloizou 1987). The central 



density of STARs is as high as ~ 100 Mq pc~ and the 
central velocity dispersion of the STARs is ~ 3 km s^^. 
This value of velocity dispersion is smaller th an the ob- 
served value (Dubath, Meylan & Mayor 1997). We need 
to follow further evolution of the star clusters for proper 
comparison, which is not feasible with present code. The 
number of STAR particles at t = 10 Myr is 2500. 

For Ri — 200 pc and 300 pc, the overall evolution is 
similar to the case for Ri = 150 pc and the stellar masses 
at t = 10 Myr are - 8.5 x 10* M© and - 5.6 x 10* Mq, 
respectively. The bound stellar mass are ^ 5 x 10* Mq for 
both case. These results are summarized in Table |[ The 
initial concentration correlates with the final stellar mass 
and the final concentration of stellar system. In order for 
the stellar system as massive as 10^ Mq to form, the PGC 
should be as compact as Ri < 200 pc (pc > 10~^^ g cm~'^) 
for the metal-free condition. 

4.2. Evolution of the metal-rich PGC 

The chemical composition of stars in globular clusters is 
one of the most crucial quantities to constrain the model 
for the globular cluster formation. The assumption that 
PGC initially has some heavy elements is quite reasonable. 
Actually, all globular clusters in our galaxy have some met- 
als of [Fe/H] = -2.25 - 0. The metallicity distribu tion of 



the g lobular clusters shows bimodal distributions (Harrii 



1991). From this fact, we choose the initial metallicity of 
the cloud ranging from [Fe/H] = —2 to 0. For comparison, 
we will evolve the lower metallicity ([Fe/H] = —3) cloud. 
The results presented in this section are obtained for case 
B. 

Using the same initial conditions presented in section 
4.1, we evolve the PGC for different metallicity. For the 
initial radius of the cloud, we set Ri = 300 pc. The first 
star formation occurs in the central region and the STARs 
heat up the surrounding matter to gradually increase the 
temperature of the central region. At i ~ 6 Myr, Type 
II supernovae begin to occur and the temperature of the 
heating region increases rapidly to T ~ 10® K. Figure || 
shows the evolutionary changes in the central tempera- 
ture (left panel) and gas density (right panel) for [Fe/H] 
= —2. The evolution for [Fe/H] = — 2 is very similar to 
the evolution of the metal-free cloud. After t ~ 6 Myr, 
the central high density region makes the shell structure 
as in the metal-free case. At t = 10 Myr, the bound stellar 
mass reaches ~ 1.0 x 10^ Mq for [Fe/H] = -2. 

For higher metallicity, i.e., [Fe/H] > — 1, details of the 
evolution are somewhat different. We compare the central 
temperature evolution for different metallicity in Figure 



10. Due to the different cooling rate, the initial decrease 
in the temperature is larger for higher metallicity. This dif- 
ference leads to different star formation history as shown 
in Figure ^ For [Fe/H] = 0, the initial star burst is 
very intense and then the SFR decrease sharply becoming 
lower than the low metallicity case after t = 1 Myr. This 
is because heating due to the stellar winds much larger 
for higher metallicity (see Section 2.3.1). The rise in the 
SFR after t = 6 Myr for [Fe/H] = is caused by Type II 
supernovae. 

The SFR for [Fe/H] = — 1 is almost constant during the 
evolution. For [Fc/H] — —2, the first star formation oc- 
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curs at i 0.5 Myr because of the lower cooling rates. The 
SFR after t = 2 Myr is almost constant but lower than for 
[Fe/H] = -1. 

The bound stellar mass at t = 10 Myr is ~ 1.5 x 10^ 
Mq and 1.0 X 10^ Mq for [Fe/H] = -1 and 0, respectively. 
The reason for such a metallicity dependence is summa- 
rized below. The bound stellar mass, the central stellar 
density, and the central velocity dispersion are summa- 
rized in Table |l]. The results depend on the metallicity as 
follows: 

1. The initial metallicity significantly affects the star 
formation history. This difference makes the subse- 
quent evolution different. 

2. For the higher metallicity, the final stellar mass at 
t = 10 Myr is larger because of the more efficient 
cooling rate. 

3. The bound stellar mass is not an increasing function 
of the initial metallicity. This is because heating due 
to the stellar winds is larger for higher metallicity. 

4. For lower metallicity ([Fe/H] = -3), only ~ 3 x 10^ 
Mq stars form. This implies that [Fe/H] > — 2 is 
necessary to form globular clusters of ~ 10^ Mq if 
Ri - 300 pc. 

4.3. Parameter dependence 

In this subsection, we describe the dependence of the 
results on various numerical parameters, e.g., the initial 
particle number (Ni) and the feedback radius (i?/). We 
use the metal-rich cloud of [Fe/H] — —1 as a reference 
model, where Ni = 5000 and Rj = 3 pc. 

First, we describe the dependence on the initial parti- 
cle number. With larger (smaller) number of initial gas 
particles, the initial masses of the gas particles are smaller 
(larger). There is no strong dependence on the mass of the 
gas particles in our star formation recipes because Eq. ( 12 ) 
dose not include the mass of the gas particles. However, 
the maximum density (/Omax) that can be represented in 
the SPH method depends on the mass of the gas particles, 
which may produce weak dependence on N^. To confirm 
this, we evolve the model with Ni ~ 10"*. The overall evo- 
lution is almost indistinguishable to the reference model. 
The stellar mass at t = 10 Myr becomes ~ 1.2 x 10^ Mq, 
which is slightly smaller than 1.3 x 10^ Mq in the reference 
model. For Ni = 2500 (a half of the reference model], we 
obtain the stellar mass of 1.5 x 10^ Mq. In Figure |2[ we 
compare the gas density profiles at t = 10 Myr for different 
Ni. The position of the shell for the model with Ni — 10** 
is different from other models because of the smaller stel- 
lar mass at the same epoch. However, all three models 
show the similar density profiles so that we conclude that 
the dependence on the initial gas particle number is week. 
Thus, Ni ~ 5000 is sufficient for the current numerical 
model. 

Secondary, we describe the dependence on the size of 
the feedback radius Rf. In our model, the cooling and 
the star formation are suppressed for the gas particles in 
heating phase. We therefore expect that the size oi Rf 
affects the star formation history. With larger Rf = 5 pc, 
the star formation rates are smaller and the final stellar 
mass (~ 0.9 x 10^ Mq) is smaller than 1.3 x 10^ Mq in 



the reference model {Rf — 5 pc). On the other hand, we 
obtain the stellar mass of ^ 1.9 x 10^ Mq with smaller 
Rf {— 2 pc). Within the reasonable range oi Rf around 
3pc, the final stellar mass changes by a factor of 0.7 - 1.5. 

4.4. Summary of the numerical results 

The results of our calculations are summarized as fol- 
lows (see also Table : 

1. In all cases, the overall evolution is similar. Initially, 
the star burst occurs in the central region. The cen- 
tral star burst is halted by the heating due to Type 
II supernovae. The heating cause the expansion of 
the central region and forms a shell structure. The 
subsequent star formation occurs in the shell. 

2. For the metal- free collapse case, the stellar mass at 
t = 10 Myr correlates with initial concentrations. To 
obtain the star cluster as massive as globular clus- 
ters {r-^ 10^ -^^0)7 the initial concentration of the 
PGC must be large enough, i.e., Ri < 200 pc. 

3. For the metal-rich collapse, the initial metallicity sig- 
nificantly affects the evolution. To obtain the star 
cluster as massive as globular clusters, the initial 
metallicity must be as larger as [Fe/H] ~ — 2. If 
the initial metallicity is low ([Fe/H] < —2), very few 
stars form. 

4. This suggests that during the initial phase of the 
galaxy formation, i.e., when the ISM contains little 
heavy element, the formation efficiency of the glob- 
ular cluster is low. 

5. The results is not strongly affected by the initial 
number of gas particles. Ni ^ 5000 is sufficient for 
the current numerical model. The dependence on 
the feedback radius is more evident. 

5. DISCUSSION 

5.1. Star formation in the shell 

In all our numerical models, a shell-like gaseous struc- 
ture is formed. The formation of the shell-like structure 
of gases has also been reported in the simulations of the 
formation of dwarf elliptical gal axies (Mori et al. 1997; 
Mori, Yoshu, fc Nomoto 1999| ). They argued that the 



difference in the density structure between normal ellip- 
ticals (de Vaucouleurs law) and dwarf ellipticals (expo- 
nential law) can be explained by the formation of such a 
shell-like structure and the star formation in the shell. The 
most crucial difference between the normal elliptical and 
the dwarf is their mass, and the less massive galaxy is more 
strongly affected by the energy inputs from stars. Similar 
argument may be applicable to our numerical models. In 
the case of the globular cluster formation, however, the 
mass is even smaller than the dwarf galaxies. Owing to 
this fact, the effect of the energy inputs from stars is more 
drastic so that the star formation in the shell is much less 
efficient than in the dwarf galaxies. Also, the stars formed 
in the shell is not gravitationally bound due to the outward 
velocity of the shell. Thus, the bound globular cluster con- 
sists of the stars formed before the shell formation. The 
stellar density of such stars is not affected by the shell 
formation and the star formation in the shell. 
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The chemical composition of stars in globular clusters is 
one of the most crucial quantities to constrain the model 
for the globular cluster formation. The stars in a glob- 
ular clu ster have almo st the same heavy elements abun- 
dances ( [Suntzeff 1993 ). This small dispersion in metallic- 
ity (crpe/HJ) suggests that the formation period of glob- 
ular clusters is so short that the stars in globular clusters 
can be regarded to form almost simultaneously as shown 
in our numerical models. 



Brown, Burkert & Truran (1991, 1995) suggested that 
the second generation stars would form in the shell and 
the self-enrichment could occur there. In our model, the 
star formation takes place in the shell after t = 5 Myr 
(see Figure However, when the stars form in the shell, 
the shell has not been enriched with newly ejected heavy 
metal as shown in the solid line in Figure 0. This line 
shows the radius of the metal-enriched region defined as 
r = (SriiTii) / (Smi) by summing over the metal-enriched 
gas particles. The metal-enriched region expands out- 
wardly but does not reach the star forming region, which 
implies that the self-enrichment dose not occur in our 
model (we note that the star formation is suppressed in 
the gas particles near the STARs in our numerical model 
as described in section 2.3.3). In the present model, no 
external medium outside the cloud is included, because of 



technical difficulties in the SPH method (see, however, Na-. 
gasawa fc Miyama (1987| ) for a possible improvement). If 
there exists external medium outside the cloud, the den- 
sity of the shell would be higher and the star formation 
history would be different; this possibility needs further 
study to confirm. 

5.2. Failed Globular clusters 

When a PGC is an initially metal-free cloud or the ini- 
tial metallicity of a PGC is low, the resulted mass and the 
central stellar density are lower than the observed mass 
and density of the globular clusters. Such "failed" globular 
clusters might be the field halo stars. Another possibility 
is a open cluster. Typical age of open clusters in our galaxy 
is r < 10^ yr. Howe ver, there al so exist such old open clus- 
ters as r > 10 Gyr (Friel 1995). The age of the most old 
open cluster is comparable to the age of globular clusters. 
The central density of the "failed" globular cluster is as 
low as the central density of typical open clusters. Thus, 
the formation processes of such old open clusters may be 
the similar to our model of globular cluster formation but 
starting from lower initial concentration. Very old open 
clusters might be the debris of "failed" globular clusters 
and there might have existed many more open clusters at 
the formation of our galaxy. 

6. CONCLUSIONS 

For the processes of globular clusters formation, only 
the qualitative scen a rios have been discus sed previously 
dFall fc Rees 1985|; |Lin fc Murray 1996|). In this pa- 



per, we present the first attempt to simulate the globular 
cluster formation with three-dimensional hydrodynamical 
method, which includes the star formation and its feed- 
back effects. We assume that, in the collapsing galaxy, 
isothermal cold clouds form through thermal condensa- 
tions and become proto-globular clouds. We obtain the 
size of proto-globular clouds by means of the linear stabil- 
ity analysis (Figure ^) and compute the evolution of the 



inner region of the PGC starting from various initial ra- 
dius Ri. The results of our calculations are summarized 
as follows: 

1. In order for the globular cluster- like system to form 
from a metal-free PGC, the initial concentration of 
the PGC must be large enough. 

2. It is required that the metallicity of a PGC is high 
enough to produce the globular cluster-like system. 
In our calculations, the required metallicity is esti- 
mated to be [Fe/H] > -2. 

3. In all cases, the shell like structure of the gas forms. 
Although the star formation occurs in the shell, the 
self-enrichment is not seen to occur. 

Based on the earlier qualitative works and the present 
quantitative results, the processes of globular clusters for- 
mation in the proto-galaxy can be understood as follows: 

1. In the collapsing proto-galaxy, the first generation 
stars of M ~ lO'^ A/q form due to the efficient cool- 
ing by hydrogen molecules. Such population HI stars 
eject the gas with heavy elements and chemically 
contaminate the proto-galaxy gases. 

2. Such first generation stars radiate dissociative pho- 
tons, and the entire proto-galaxy is settled down to 
ionization equilibrium. 

3. With equilibrium cooling, the density perturbation 
grows due to thermal instability, and forms an 
isothermal cloud with a density structure of cx r~^. 
Such clouds are the proto-cloud of globular clusters. 

4. When the density of the PGC becomes high enough, 
the burst of star formation occurs. Some high den- 
sity clouds produce the globular clusters, and others 
may produce the field stars and/or the halo open 
clusters. 

5. During the formation of the galaxy, the formation ef- 
ficiency of the globular cluster become significantly 
large when the metallicity of the PGC become as 
large as [Fe/H] 2. 

As mentioned in section 5.1, we do not take account of 
the effect of the external medium in the present calcula- 
tions. If we include the effect of the external medium, the 
structure of the shell may be dif ferent and further star for- 
matio n may occur in that shell. Brown, Burkert & Truran 
|(1995| ) suggested that only a few supernovae per Myr is 
sufficient to reverse the contraction of the 10® Mq cloud. 
Such an effect of the external medium needs to be studied. 

The globular cluster formation is considered to oc- 
cur not only in the proto-galaxies but also in the 
present galaxies, e.g., LMC, SMC, and interacting galaxy 
NGC4038/NGC4039. In such environment, different for- 
mation processes wo uld take place, which correspond t o 
the tertiary models ( pK^umai, Basu, fc Fujimoto 1993 ). 
In the case of NGC4038/NGC4039 (the Antennae galax- 
ies), many globular clusters are being produced by galaxy 
merging. Detail formation processes of globular clusters in 
such systems needs to be studied. Also, the effect of the 
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formation of many globular clusters on the galaxy merging 
process needs a further study. 
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Table 1 

The summary of the numerical results 



r> a 
-fti 


[Fe/H]b 


Mass'^ 






150 pc 


no metal 


1.0 X 10^ Mq 


lO^-S Mq pC-3 


2.98 km s"^ 


200 pc 


no metal 


5.1 X 10" Mq 


IQi-i Mq pc-3 


1.87 km 


300 pc 


no metal 


5.5 X 104 Mq 


101-4 Mq pC-3 


3.66 km s~i 


300 pc 





1.0 X 10-^ Mq 


10^-" Mq pC-3 


3.21 km 


300 pc 


-1 


1.5 X 10^ Mq 


10^-0 Mq pc-3 


2.33 km 


300 pc 


-2 


1.0 X 10^ Mq 


lO^-^ Mq pC-3 


2.68 km s"^ 


300 pc^ 


-3 


2.9 X 10^ Mq 







Note.— ^The initial radius of the PGC. 

I'Thc initial mctallicity ([Fc/H]) of the PGC. 

'^The bound stellar mass at t = 10 Myr. 

•^The central stellar density at t = 10 Myr. 

*^The central stellar velocity dispersion at t = 10 Myr. 

^ Pc and CTc is omitted of this case. 



Table 2 

The dependence on numerical parameters 



A^i* Rf^ Mass" 



5000 3 pc 1.5 X 10^ Mq 

2500 3 pc 1.5 X 10^ Mq 

10* 3 pc 1.2 X 10^ Mq 

5000 2 pc 1.9 X 10^ Mq 

5000 5 pc 8.6 X 10* Mq 



Note. — ''Used initial gas particle number. 
I' The feedback radius. 
'^The stellar mass at i = 10 Myr. 
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log T (K) 



Fig. 1. The various cooling rates used in SPH code. The solid and dotted show the cooling rates for cases A and C, respectively. The 
dashed lines show the cooling rates for case B. From the top to the bottom line, each line corresponds to [Fe/H] = —1, —2 and —3, respectively. 
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log T (K) 

Fig. 2. — The ratio between the recombination time of hydrogen and the coohng time (trccom/tcool) a function of the temperature. 




5x10^ 10^ 1.5x10^ 

log T (K) 

Fig. 3. — The estimated scale length of a proto-globular cloud corresponding to the maximum growth rate in Eq. (p^). The solid, dashed, 
and dotted lines correspond to po = 1.0 X 10"'^^, 1.0 X 10"^*, and 1.0 X lO^^'' g cm~^, respectively. 
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Fig. 4. — The initial gas density (g cm ^) profile for the King profile case. The solid, dashed, and dotted lines show the profiles for Ri 
150, 200, and 300 pc, respectively. 
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Fic. 5. The evolutionary changes in the central temperature (left panel) and gas density (right panel) are shown for the metal-free 
collapse with Ri = 150 pc. 
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Fig. 6. — The gas density profile as a function of radius (pc). left: The soUd, dashed and dotted hnes respectively shows the profile at t = 
0, 2, 4 Myr. right: The solid lines with the decreasing central density correspond to t = 5, 6, 7, 8 Myr. 
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Fig. 7. — The radius (pc), at which stars form, as a function of time (yr) for the metal-free collapse with i?j = 150 pc. 
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Fig. 8. — The projected position of particles at t = 10 Myr for metal-free collapse with Ri = 150 pc (left panel). The stellar density profile 
at t = 10 Myr for the same model. 
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Fig. 9. — The evolutionary changes in the central temperature (left panel) and gas density (right panel) are shown for the collapse of 
[Fe/H] = -2 sphere. 
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Fig. 10. — The evolutionary changes in the central temperature for the collapse of [Fe/H] 
line) sphere. 
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Fig. 11. — The comparison of SFR as a function of time for [Fe/H] = —2 (solid line), —1 (dashed line) and (dotted line). 
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Fig. 12. — The comparison between the reference model (sohd hne) and other model. The dashed line corresponds to the model with iVj ~ 
2500. The dashed line corresponds to the model with ATj ~ 10000. 



